
#include "skinning.h"

//~ namespace MLS {

class compareuv {
public:
   bool operator() (const MeshWithSkeleton::Vertex_handle &u, 
                    const MeshWithSkeleton::Vertex_handle &v) const {
      return u->d > v->d;
   }
};

Skinning::Skinning(MeshWithSkeleton *m) {
   
   model = m;
   stack_search.resize(m->nBones()); // 
   lflist_bones.resize(model->nJoints());
   
   init_vars();
}

Skinning::Skinning() {
   
   model = 0;
   init_vars();
   //~ lflist_bones.resize(model->nJoints());
}

void Skinning::init_vars() {
   
   diffusion_function_option=0;
   diffusion_amplitude=0.1;
   diffusion_precision=100;
   normalized_gaussian_table.resize(diffusion_precision);
   standard_deviation = diffusion_amplitude/3.0;
   
   limit_weights=20;
   list_size=100000;
   
   //~ inter_test_limit = 1.0; 
   
   //~ stack_search.resize(m->nBones()); // 
   
   num_iterations = 20;
   num_expansions = 20;
}

void Skinning::execute(MeshWithSkeleton *m) {
   
   cout<<"execute():list_size: "<<list_size<<endl;
   
   assert(m != 0);
   model = m;
   stack_search.resize(m->nBones()); // 
   
   execute();
}

void Skinning::execute() { // procedure skinning
   
   assert(model != 0);
   
   cout<<"skinning begin!"<<endl;
   cout<<" num_iterations: "<<num_iterations<<endl;
   cout<<" num_expansions: "<<num_expansions<<endl;
   cout<<" diffusion_amplitude: "<<diffusion_amplitude<<endl;
   
   //~ // Reset variables for skinning
   stack_search.clear();
   stack_search.resize(model->nBones());
   //~ component_list.clear();
   lflist_bones.clear();
   lflist_bones.resize(model->nBones()); 
   
   cout<<"xxx model->nBones(): "<<model->nBones()<<endl;
   
   
   MeshWithSkeleton::Vertex_iterator vertex = model->vertices_begin();
   /// Initialize the joints list for each vertex
   MeshWithSkeleton::Vertex_iterator vi = model->vertices_begin();
   
   while (vi != model->vertices_end()) {
      assert(vi->rigging_group != -1);
      // clear old data
      //~ vi->reset();
      
      //~ vi->assigned = false;
      //~ vi->jointsweights.clear();
      //~ vi->jointsweights.resize(model->nJoints(), 0.0);
      //~ vi->jointsweights[vi->rigging_group] = 1.0;
      
      vi->weightlf[0] = 0;
      vi->weightlf[1] = 0;
      vi->assigned = false;
      vi->bonesweights.clear();
      vi->bonesweights.resize(model->nBones(), 0.0);
      vi->bonesweights[vi->rigging_group] = 1.0;
      
      ++vi;
   }
   
   
   
   
   /// ATTTENTION ATTENTIONNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
   model->computeLaplacianParameters();
   
   
   cout<<" find connected component ... begin ..."<<endl;
   /// skinning 1
   //~ vertex = model->vertices_begin();
   //~ int current_group=0, current_component=0;
   //~ for (int i=0; i<model->nvertices(); ++i) {
      //~ if (vertex[i].assigned) continue;
      //~ find_connected_component(i, current_component);
      //~ ++current_component;
   //~ }
   
   vertex = model->vertices_begin();
   for (int i=0; i<model->nvertices(); ++i) {
      if (vertex[i].assigned) continue;
      find_connected_component_lf(i, vertex[i].rigging_group);
   }
   cout<<" find connected component ... done!"<<endl;
   
   cout<<" perform propagation ... begin ..."<<endl;
   /// skinning 1
   //~ vertex = model->vertices_begin();
   //~ current_group = current_component;
   //~ for (int i=0; i<model->nJoints(); ++i) {
      //~ propagate_rigging_group(i, current_group);
      //~ ++current_group;
   //~ }
   
   /// skinning nuevo
   int current_group = model->nBones();
   for (int i=0; i<model->nBones(); ++i) {
      //~ cout<<"wld para bone!: "<<i<<endl;
      weightsLaplacianDiffusion(i, current_group);
      ++current_group;
   }
   cout<<" perform propagation ... done!"<<endl;
   
   /// Weights normalization
   /// skinning 1
   //~ cout<<" weights normalization... begin ..."<<endl;
   //~ for (int i=0; i<model->nvertices(); ++i) {
      //~ Number sumtmp = 0;
      //~ assert(vertex[i].joints.size() != 0);
      //~ for (int j=0; j<vertex[i].joints.size(); ++j) {
         //~ sumtmp += vertex[i].jointsweights[j];
      //~ }
      //~ for (int j=0; j<vertex[i].joints.size(); ++j) {
         //~ vertex[i].jointsweights[j]/=sumtmp;
      //~ }
   //~ }
   
   cout<<" analise dos pesos!"<<endl;
   
   vertex = model->vertices_begin();
   for (int i=0; i<model->nvertices(); ++i) {
      Number sumtmp = 0;
      //~ assert(vertex[i].joints.size() != 0);
      //~ if (vertex[i].joints.size() == 0) {
         //~ vertex[i].joints.push_back(vertex[i].rigging_group);
         //~ vertex[i].jointsweights.push_back(1.0);
      //~ }
      
      for (int j=0; j<vertex[i].bonesweights.size(); ++j) {
         sumtmp += vertex[i].bonesweights[j];
      }
      //~ if (sumtmp > 1.01) {
         //~ cout<<"i: "<<i<<" "<<sumtmp<<" size: "<<vertex[i].bonesweights.size()<<endl;
         //~ for (int j=0; j<vertex[i].bonesweights.size(); ++j) {
            //~ cout<<"~:: "<<vertex[i].bonesweights[j]<<endl;
         //~ }
      //~ }
      //~ if (sumtmp < 0.99) {
         cout<<"i: "<<i<<" "<<sumtmp<<endl;
         for (int j=0; j<vertex[i].bonesweights.size(); ++j) {
            cout<<"~:: "<<vertex[i].bonesweights[j]<<endl;
         }
      //~ }
      
      
      //~ for (int j=0; j<vertex[i].joints.size(); ++j) {
         //~ sumtmp += vertex[i].jointsweights[j];
      //~ }
      //~ if (sumtmp > 1.01) {
         //~ cout<<"i: "<<i<<" "<<sumtmp<<" size: "<<vertex[i].joints.size()<<endl;
         //~ for (int j=0; j<vertex[i].joints.size(); ++j) {
            //~ cout<<"~:: "<<vertex[i].joints[j]<<" :: "<<vertex[i].jointsweights[j]<<endl;
         //~ }
      //~ }
      //~ if (sumtmp < 0.99) {
         //~ cout<<"i: "<<i<<" "<<sumtmp<<" size: "<<vertex[i].joints.size()<<endl;
         //~ for (int j=0; j<vertex[i].joints.size(); ++j) {
            //~ cout<<"~:: "<<vertex[i].joints[j]<<" :: "<<vertex[i].jointsweights[j]<<endl;
         //~ }
      //~ }
   }
   cout<<" weights normalization... done!"<<endl;
   
   
   //~ vertex = model->vertices_begin();
   //~ for (int i=0; i<model->nvertices(); ++i) {
      //~ cout<<"i: "<<i<<" ::: \n ";
      //~ for (int j=0; j<vertex[i].joints.size(); ++j) {
         //~ cout<<"   vertex[i].joints: "<<vertex[i].joints[j]<<" :: "<<vertex[i].jointsweights[j]<<endl;
      //~ }
      
      //~ cout<<"<> "<<i<<" <> "<<vertex[i].rigging_group<<endl;
   //~ }
   
   
}

/// Computes the normalized_gaussian_table
void Skinning::compute_normalized_gaussian_table() {
   
   Number temp = diffusion_amplitude/(diffusion_precision*standard_deviation);
   Number step2 = temp*temp;
   temp = diffusion_amplitude/standard_deviation;
   Number min_exp = exp(-(temp*temp));
   for (int i=0; i<diffusion_precision; ++i)
      normalized_gaussian_table[i] =(exp( -(i*i*step2)) - min_exp)/(1.0-min_exp);
}

Number Skinning::diffusion_function(int diffusion_option, Number dist) {
   
   switch (diffusion_option) {
      case 0: {
         Number sq_root = 1.0 - dist/diffusion_amplitude;
         return sq_root*sq_root;
      }
      case 1: {
         return normalized_gaussian_table[(int)(dist/diffusion_amplitude*diffusion_precision)];
      }
   }
}

Number Skinning::length(int id_he) const {
   
   MeshWithSkeleton::Halfedge_iterator edge = model->halfedges_begin();
   const Point3 &pA = edge[id_he].vertex()->point();
   const Point3 &pB = edge[id_he].opposite()->vertex()->point();
   const Number &vx = pB[0]-pA[0];
   const Number &vy = pB[1]-pA[1];
   const Number &vz = pB[2]-pA[2];
   
   return sqrt(vx*vx + vy*vy + vz*vz);
}

void Skinning::find_connected_component_lf(int v, int current_component) {
   
   //~ cout<<"find_connected_component_lf: "<<v<<" cc: "<<current_component<<endl;
   
   int current_edge, current_vertex, edge_ant, edge_ini;
   Number distance_to_bone;
   struct edge_pair {
      int e_ini;
      int e_term;
   };
   edge_pair list[list_size];
   
   MeshWithSkeleton::Vertex_iterator   vertex = model->vertices_begin();
   MeshWithSkeleton::Halfedge_iterator edge = model->halfedges_begin();
   
   list[0].e_ini = vertex[v].halfedge()->hedgeid; // 
   list[0].e_term = vertex[v].halfedge()->hedgeid;
   int current_bone = vertex[v].rigging_group;
   int n_listed = 1;
   int n_explored = 0;
   
   //~ cout<<"::rigging_group: "<<vertex[v].rigging_group<<endl;
   
   while (n_listed > n_explored) {
      //~ if (current_bone==1)
      //~ cout<<" ..while: "<<n_listed<<" :: "<<n_explored<<endl;
      current_edge = list[n_explored].e_ini;
      edge_ant = edge[current_edge].next()->hedgeid;
      do {
         current_vertex = edge[current_edge].vertex()->id();
         //~ cout<<" .... do: "<<current_vertex<<endl;
         //~ cout<<" .... ##: "<<vertex[current_vertex].rigging_group<<" "<<current_bone<<endl;
         if (vertex[current_vertex].rigging_group == current_bone) {
            //~ cout<<" ..ifififif: "<<endl;
            if ((vertex[current_vertex].assigned)) { 
               current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
               edge_ant = edge[current_edge].next()->hedgeid;
            }
            else {
               vertex[current_vertex].assigned = true;
               vertex[current_vertex].weightlf[0] = 1;
               vertex[current_vertex].weightlf[1] = 1;
               list[n_listed].e_ini = edge[edge[edge_ant].opposite()->hedgeid].next()->hedgeid;
               current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
               edge_ant = edge[current_edge].next()->hedgeid;
               list[n_listed].e_term = edge[edge_ant].opposite()->hedgeid;
               n_listed = (n_listed%list_size) + 1;
            }
         }
         else {
            //~ cout<<" >>>>..else: "<<current_vertex<<" "<<current_edge<<endl;
            if (vertex[current_vertex].curr_comp != current_component) {
               //~ cout<<" #### if 111111111: "<<endl;
               vertex[current_vertex].curr_comp = current_component;
               //~ cout<<" #### if 222222: "<<endl;
               current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
               //~ cout<<" #### if 333333333: "<<endl;
               edge_ant = edge[current_edge].next()->hedgeid;
                   
               //~ cout<<" #### if 444444: "<<endl;
               
                   
               lflist_element record(current_vertex, vertex[current_vertex].halfedge()->opposite()->hedgeid, 
                                     vertex[current_vertex].halfedge()->opposite()->hedgeid);
               //~ cout<<" #### if 55555555555: "<<vertex[v].rigging_group<<endl;
               lflist_bones[vertex[v].rigging_group].push_back(record);
               //~ cout<<" #### if 66666: "<<endl;
            }
            else {
               //~ cout<<" #### else "<<endl;
               current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
               edge_ant = edge[current_edge].next()->hedgeid;
            }
         }
      }
      while(current_edge != list[n_explored].e_term);
      n_explored = (n_explored%list_size) + 1;
   }
   
}

void Skinning::weightsLaplacianDiffusion(int bone, int current_group) {
   
   //~ cout<<"lflist_bones[bone].size(): "<<lflist_bones[bone].size()<<endl;
   //~ cout<<"   weightsLaplacianDiffusion()"<<endl;
   
   int nlisted = lflist_bones[bone].size();
   int nlistedoriginal = nlisted;
   int nlistedprevious = 0;
   int neww = 1, oldw = 0; // new and old weights
   float LFSTOPDELTA = 0; /// /////////////////////////////////////////
   
   Number min_weight;
   int i_min_weight; 
   lflist_element record;
   int current_edge, edge_ant, current_vertex;
   
   MeshWithSkeleton::Vertex_iterator   vertex = model->vertices_begin();
   MeshWithSkeleton::Halfedge_iterator edge = model->halfedges_begin();
   
   //~ cout<<"   1111111 nlisted ant: "<<nlisted<<" bone: "<<bone<<" current_group: "<<current_group<<endl;
   
   for (int j=0; j<nlisted; ++j) {
      record = lflist_bones[bone][j];
      int initial_vertex = record.vertex;
      while (vertex[initial_vertex].curr_comp == current_group) {
         //~ cout<<"zzzzzzz: "<<initial_vertex<<" "<<current_group<<" j: "<<j<<" "<<nlisted<<" :: "<<lflist_bones[bone].size()<<endl;
         lflist_bones[bone][j] = lflist_bones[bone][nlisted-1];
         //~ cout<<":> "<<lflist_bones[bone].size()<<endl;
         --nlisted;
         if (j==nlisted) break;
         record = lflist_bones[bone][j];
         initial_vertex = record.vertex;
      }
      if (j==nlisted) break;
      
      vertex[initial_vertex].curr_comp = current_group;
      vertex[initial_vertex].weightlf[0] = 0;
      vertex[initial_vertex].weightlf[1] = 0;
      current_edge = record.e_init;
      edge_ant = edge[current_edge].next()->hedgeid;
      do {
         //~ cout<<"current_vertex: "<<current_vertex
         //~ <<" .curr_comp: "<<vertex[current_vertex].curr_comp
         //~ <<" rigging_group "<<vertex[current_vertex].rigging_group<<endl;
         current_vertex = edge[current_edge].vertex()->id();
         if ((vertex[current_vertex].curr_comp != current_group) and 
            (vertex[current_vertex].rigging_group == bone)) {
            vertex[current_vertex].weightlf[0] = 1;
            vertex[current_vertex].weightlf[1] = 1;
            vertex[current_vertex].curr_comp = current_group;
               
            //~ cout<<"current_vertex: "<<current_vertex<<" current_group: "<<current_group<<" rigging_group "<<vertex[current_vertex].rigging_group<<endl;
            
            //~ int edge_ini = edge[edge[edge_ant].opposite()->hedgeid].next()->hedgeid;
            current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
            edge_ant = edge[current_edge].next()->hedgeid;
            //~ int edge_term = edge[edge_ant].opposite()->hedgeid;
            //~ lflist_element record(current_vertex, edge_ini, edge_term);
               
               
            lflist_element newrecord(current_vertex, vertex[current_vertex].halfedge()->opposite()->hedgeid, 
               vertex[current_vertex].halfedge()->opposite()->hedgeid);
            lflist_bones[bone].push_back(newrecord);
         }
         else {
            current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
            edge_ant = edge[current_edge].next()->hedgeid;
         }
      }
      while (current_edge != record.e_term);
   }
   
   for (int k=nlisted; k<nlistedoriginal; ++k) {
      lflist_bones[bone][k] = lflist_bones[bone][lflist_bones[bone].size()-1];
      lflist_bones[bone].pop_back();
   }
   
   nlisted = lflist_bones[bone].size();
   nlistedoriginal = nlisted;
   
   //~ cout<<"   nlisted dep: "<<nlisted<<endl;
   //~ cout<<"   22222222 nlisted ant: "<<nlisted<<endl;
   
   int i=0; bool stop_flag = true;
   for (i=1; i<num_expansions; ++i) {
      stop_flag = true;
      
      for (int j=nlistedprevious; j<nlisted; ++j) {
         record = lflist_bones[bone][j];
         current_edge = record.e_init; 
         edge_ant = edge[current_edge].next()->hedgeid;
         do { 
            current_vertex = edge[current_edge].vertex()->id();
            if (vertex[current_vertex].curr_comp != current_group) {
               vertex[current_vertex].curr_comp = current_group;
               
               int edge_ini = edge[edge[edge_ant].opposite()->hedgeid].next()->hedgeid;
               current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
               edge_ant = edge[current_edge].next()->hedgeid;
               int edge_term = edge[edge_ant].opposite()->hedgeid;
               if (vertex[current_vertex].rigging_group == bone) {
                  vertex[current_vertex].weightlf[oldw] = 1;
                  vertex[current_vertex].weightlf[neww] = 1;
               }
               else {
                  vertex[current_vertex].weightlf[oldw] = 0;
                  vertex[current_vertex].weightlf[neww] = 0;
               }
               
               
               lflist_element newrecord(current_vertex, edge_ini, edge_term);
               lflist_bones[bone].push_back(newrecord);
            }
            else {
               current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
               edge_ant = edge[current_edge].next()->hedgeid;
            }
         }
         while (current_edge != record.e_term);
         
         int idv = lflist_bones[bone][j].vertex;
         computeLaplacianFilter(vertex[idv].halfedge()->vertex(), bone, neww);
      }
      
      //~ cout<<"i before : "<<i<<endl;
      
      if (nlistedprevious == 0) {
         neww = oldw;
         oldw = 1-oldw;
         nlistedprevious = nlisted;
         nlisted = lflist_bones[bone].size();
         continue;
      }
      
      for (int j=0; j<nlistedprevious; ++j) {
         int idv = lflist_bones[bone][j].vertex;
         computeLaplacianFilter(vertex[idv].halfedge()->vertex(), bone, neww);
         if (stop_flag and ((abs(vertex[idv].weightlf[neww]-vertex[idv].weightlf[oldw]) > LFSTOPDELTA))) {
            stop_flag = false;
         }
      }
      //~ for (int j=nlistedoriginal; j<nlistedprevious; ++j) {
         //~ int idv = lflist_bones[bone][j].vertex;
         //~ computeLaplacianFilter(vertex[idv].halfedge()->vertex(), bone, neww);
      //~ }
      
      if (stop_flag) break;
      neww = oldw;
      oldw = 1-oldw;
      nlistedprevious = nlisted;
      nlisted = lflist_bones[bone].size();
   }
   
   //~ cout<<"   33333333 nlisted ant: "<<nlisted<<endl;
   
   /// ////////
   if (not stop_flag) {
      for (int i=num_expansions; i<num_iterations; ++i) {
         for (int j=0; j<nlistedprevious; ++j) {
            int idv = lflist_bones[bone][j].vertex;
            computeLaplacianFilter(vertex[idv].halfedge()->vertex(), bone, neww);
            if (stop_flag and ((abs(vertex[idv].weightlf[neww]-vertex[idv].weightlf[oldw]) > LFSTOPDELTA))) {
               stop_flag = false;
            }
         }
         //~ for (int j=nlistedoriginal; j<nlistedprevious; ++j) {
            //~ int idv = lflist_bones[bone][j].vertex;
            //~ computeLaplacianFilter(vertex[idv].halfedge()->vertex(), bone, neww);
         //~ }
         
         if (stop_flag) break;
            
         neww = oldw;
         oldw = 1-oldw;
      }
   }
   
   //~ cout<<"bone: "<<bone<<endl;
   //~ cout<<"   44444444 nlisted ant: "<<nlisted<<endl;
   
   for (int i=0; i<nlistedprevious; ++i) {
      int idv = lflist_bones[bone][i].vertex;
      //~ cout<<"~vertex[idv].weightlf[0]: "<<vertex[idv].weightlf[0]<<endl;
      //~ cout<<"~~vertex[idv].weightlf[1]: "<<vertex[idv].weightlf[1]<<endl;
      Number weight_tmp = 0.5*(vertex[idv].weightlf[0] + vertex[idv].weightlf[1]);
      //~ cout<<"weight_tmp: "<<weight_tmp<<endl;
      if (weight_tmp < 1e-3) continue;
      //~ vertex[current_vertex].weightlf[0] = 0;
      //~ vertex[current_vertex].weightlf[1] = 0;
      
      vertex[idv].bonesweights[bone] = weight_tmp;
      
      /*
      if (vertex[idv].joints.size() >= limit_weights) {
         min_weight = 1E8;
         for (int j=0; j<vertex[idv].joints.size(); ++j) {
            if (vertex[idv].jointsweights[j] < min_weight) {
               min_weight = vertex[idv].jointsweights[j];
               i_min_weight = j;
            }
         }
         if (weight_tmp > min_weight) {
            //~ vertex[idv].jointsweights[i_min_weight] = weight_tmp;
            //~ vertex[idv].joints[i_min_weight] = bone;
            
            vertex[idv].bonesweights[i_min_weight] = weight_tmp;
         }
      }
      else {
         //~ cout<<"weight_tmp: "<<weight_tmp<<endl;
         //~ vertex[idv].jointsweights.push_back(weight_tmp);
         //~ vertex[idv].joints.push_back(bone);
         
         //~ cout<<"->vertex[record.vertex].joints.size(): "<<vertex[record.vertex].joints.size()<<endl;
      }
      */
   }
}

void Skinning::computeLaplacianFilter(Vertex_handle v, int bone, int s) {
   
   //~ cout<<"LFilter: "<<v->id()/3<<" "<<s<<" "<<v->neighbors.size()<<endl;
   int t = 1-s;
   Number wtmp = 0.0;
   for (int i=0; i<v->neighbors.size(); ++i) {
      wtmp += v->cots[i]*(v->neighbors[i]->weightlf[t]);
      //~ cout<<"v->cots[i]: "<<v->cots[i]<<endl;
      //~ cout<<"v->neighbors[i]->weightlf[t]: "<<v->neighbors[i]->weightlf[t]<<endl;
      //~ cout<<"wtmp: "<<wtmp<<endl;
   }
   
   //~ cout<<"v->lfdenom: "<<v->lfdenom<<endl;
   //~ cout<<"v->lfkterm: "<<v->lfkterm<<endl;
   
   
   if (v->rigging_group == bone) {
      v->weightlf[s] = wtmp*v->lfdenom + v->lfkterm;
   }
   else {
      v->weightlf[s] = wtmp*v->lfdenom;
   }
   //~ cout<<":: v->weightlf[s]: "<<v->weightlf[s]<<" id: "<<v->id()/3<<endl;
   
}

void Skinning::find_connected_component(int v, int current_component) {
   
   //~ int current_edge, current_vertex, edge_ant, edge_ini;
   //~ Number distance_to_bone;
   //~ struct edge_pair {
      //~ int e_ini;
      //~ int e_term;
   //~ };
   //~ edge_pair list[list_size];
   
   //~ MeshWithSkeleton::Vertex_iterator   vertex = model->vertices_begin();
   //~ MeshWithSkeleton::Halfedge_iterator edge = model->halfedges_begin();
   
   //~ list[0].e_ini = vertex[v].halfedge()->hedgeid; // 
   //~ list[0].e_term = vertex[v].halfedge()->hedgeid;
   //~ int current_bone = vertex[v].rigging_group;
   //~ int n_listed = 1;
   //~ int n_explored = 0;
   
   //~ cout<<"::rigging_group: "<<vertex[v].rigging_group<<endl;
   
   //~ while (n_listed > n_explored) {
      //~ if (current_bone==1)
         //~ cout<<" .... while"<<endl;
      //~ current_edge = list[n_explored].e_ini;
      //~ edge_ant = edge[current_edge].next()->hedgeid;
      //~ do {
         //~ current_vertex = edge[current_edge].vertex()->id();
         //~ if (vertex[current_vertex].rigging_group == current_bone) {
            //~ if ((vertex[current_vertex].assigned)) { 
               //~ current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
               //~ edge_ant = edge[current_edge].next()->hedgeid;
            //~ }
            //~ else {
               //~ vertex[current_vertex].assigned = true;
               //~ list[n_listed].e_ini = edge[edge[edge_ant].opposite()->hedgeid].next()->hedgeid;
               //~ current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
               //~ edge_ant = edge[current_edge].next()->hedgeid;
               //~ list[n_listed].e_term = edge[edge_ant].opposite()->hedgeid;
               //~ n_listed = (n_listed%list_size) + 1;
            //~ }
         //~ }
         //~ else {
            //~ Number distance_to_component = length(current_edge);
            //~ if (current_bone==1) {
               //~ cout<<" .>. "<<distance_to_component<<" "<<vertex[current_vertex].dist_comp<<endl;
               //~ cout<<vertex[current_vertex].curr_comp<<" "<<current_component<<endl;
            //~ }
            //~ if ((distance_to_component < diffusion_amplitude)  and
                //~ ((vertex[current_vertex].curr_comp != current_component) or 
                //~ (vertex[current_vertex].dist_comp > distance_to_component))) {
                   
               //~ if (current_bone==1)
                  //~ cout<<"else if"<<endl;
                   
               //~ vertex[current_vertex].curr_comp = current_component;
               //~ vertex[current_vertex].dist_comp = distance_to_component;
               //~ edge_ini = edge[edge[edge_ant].opposite()->hedgeid].next()->hedgeid;
                   
               //~ current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
               //~ edge_ant = edge[current_edge].next()->hedgeid;
                   
               //~ int edge_term = edge[edge_ant].opposite()->hedgeid;
               //~ stack_element record(current_vertex, edge_ini, edge_term, distance_to_component);
               //~ cout<<"vertex[v].rigging_group: "<<vertex[v].rigging_group<<endl;
                  //~ (record);
               //~ stack_search[vertex[v].rigging_group].push(record);
            //~ }
            //~ else {
               //~ current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
               //~ edge_ant = edge[current_edge].next()->hedgeid;
            //~ }
         //~ }
      //~ }
      //~ while(current_edge != list[n_explored].e_term);
      //~ n_explored = (n_explored%list_size) + 1;
   //~ }
}

void Skinning::propagate_rigging_group(int bone, int current_group) {
   /*
   MeshWithSkeleton::Vertex_iterator   vertex = model->vertices_begin();
   MeshWithSkeleton::Halfedge_iterator edge = model->halfedges_begin();
   
   Number min_weight, eval_func, distance_to_group;
   int i_min_weight; 
   stack_element record;
   int current_edge, edge_ant, current_vertex;
   
   while (!(stack_search[bone].empty())) {
      record = stack_search[bone].top();
      stack_search[bone].pop();
      
      if (vertex[record.vertex].curr_comp != current_group) {
         vertex[record.vertex].curr_comp = current_group;
         vertex[record.vertex].dist_comp = record.dist;
      }
      else {
         if (vertex[record.vertex].dist_comp < record.dist) { 
            continue;
         }
         else {
            vertex[record.vertex].dist_comp = record.dist;
         }
      }
      
      eval_func = diffusion_function(diffusion_function_option, record.dist);
      if (vertex[record.vertex].joints.size() >= limit_weights) {
         min_weight = 1E8;
         for (int i=0; i<vertex[record.vertex].joints.size(); ++i) {
            if (vertex[record.vertex].jointsweights[i] < min_weight) {
               min_weight = vertex[record.vertex].jointsweights[i];
               i_min_weight = i;
            }
         }
         if (eval_func > min_weight) {
            vertex[record.vertex].jointsweights[i_min_weight] = eval_func;
            vertex[record.vertex].joints[i_min_weight] = bone;
         }
      }
      else {
         if (vertex[record.vertex].joints.back() == bone ) {
         }
         vertex[record.vertex].jointsweights.push_back(eval_func);
         vertex[record.vertex].joints.push_back(bone);
      }
      
      current_edge = record.e_init; 
      edge_ant = edge[current_edge].next()->hedgeid;
      do { 
         current_vertex = edge[current_edge].vertex()->id()/3;
         Number distance_to_group = record.dist + length(current_edge);
         if((distance_to_group < diffusion_amplitude) and 
            ((vertex[current_vertex].curr_comp != current_group) or 
            (vertex[current_vertex].dist_comp > distance_to_group)) and 
            (vertex[current_vertex].rigging_group != bone)) {
            
            vertex[current_vertex].curr_comp = current_group;
            vertex[current_vertex].dist_comp = distance_to_group;
            
            int edge_ini = edge[edge[edge_ant].opposite()->hedgeid].next()->hedgeid;
            current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
            edge_ant = edge[current_edge].next()->hedgeid;
            int edge_term = edge[edge_ant].opposite()->hedgeid;
            stack_element record(current_vertex, edge_ini, edge_term, distance_to_group);
            stack_search[bone].push(record);
         }
         else {
            current_edge = edge[edge[current_edge].opposite()->hedgeid].next()->hedgeid;
            edge_ant = edge[current_edge].next()->hedgeid;
         }
      }
      while (current_edge != record.e_term);
   }
   */
}

//~ } // end MLS namespace
